Imaging of scattering media using relative detector values

ABSTRACT

A method and system for imaging of a scattering target medium using a modified perturbation formulation of a radiation transport equation where normalized measured values are used to recover a relative difference in absorption and/or scattering properties based on the normalized measured values with respect to a reference medium. The modified perturbation formulation provides enhanced stability, reduces the sensitivity of solution to variations between the target and reference media, produces solutions having physical units and reduces the need for absolute detector calibration. Moreover, the modified perturbation equation lends itself to the detection and imaging of dynamic properties of the scattering medium.

This application claims the benefit under 35 U.S.C. §120 of prior U.S.Provisional Patent Application Ser. Nos. 60/153,769 filed Sep. 14, 1999,entitled TOMOGRAPHY IN A SCATTERING MEDIUM, 60/153,926 filed Sep. 14,1999, entitled DYNAMIC TOMOGRAPHY IN A SCATTERING MEDIUM and 60/154,099filed Sep. 15, 1999, entitled DYNAMIC TOMOGRAPHY IN A SCATTERING MEDIUM.

This application is related to copending application Ser. No.10/088,254, issued as U.S. Pat. No. 6,795,195 on Sep. 21, 2004, filed onthe same date as this application, entitled “SYSTEM AND METHOD FORTOMOGRAPHIC IMAGING OF DYNAMIC PROPERTIES OF A SCATTERING MEDIUM” byinventors R. Barbour and C. Schmitz and is hereby incorporated byreference (hereinafter the “Barbour 4147PC1 application”).

This application is also related to copending application Ser. No.10/088,190, filed on the same date as this application, entitled “METHODAND SYSTEM FOR IMAGING THE DYNAMICS OF A SCATTERING MEDIUM” by inventorR. Barbour and is hereby incorporated by reference (hereinafter the“Barbour 4147PC2 application”).

This invention was made with U.S. Government support under contractnumber CA-RO166184-02A, awarded by the National Cancer Institute. TheU.S. Government has certain rights in the invention.

FIELD OF THE INVENTION

The invention relates generally to imaging in a scattering medium, andmore particularly, to a method using a novel modification to theperturbation formulation of the radiation transport inverse problem todetermine relative changes in the absorption and/or scatteringproperties of the medium based on relative changes in measured energy.

BACKGROUND OF THE INVENTION

Many techniques and systems have been developed to image the interiorstructure of a turbid medium through the measurement of energy thatbecomes scattered upon being introduced into a medium. Typically, asystem for imaging based on scattered energy detection includes a sourcefor directing energy into a target medium and a plurality of detectorsfor measuring the intensity of the scattered energy exiting the targetmedium at various locations with respect to the source. Based on themeasured intensity of the energy exiting the target medium, it ispossible to reconstruct an image representing the cross-sectionalscattering and/or absorption properties of the target. Exemplary methodsand systems are disclosed in Barbour et al., U.S. Pat. No. 5,137,355,entitled “Method of Imaging a Random Medium,” (hereinafter the “Barbour'355 patent”), Barbour, U.S. Pat. No. 6,081,322, entitled “NIR ClinicalOpti-Scan System,” (hereinafter the “Barbour '322 patent”), the Barbour4147PC1 application, and the Barbour 4147PC2 application.

Imaging techniques based on the detection of scattered energy arecapable of measuring the internal absorption, scattering and otherproperties of a medium using sources whose penetrating energy is highlyscattered by the medium. Accordingly, these techniques permit the use ofwavelengths and types of energy not suitable for familiar transmissionimaging techniques. Thus they have great potential for detectingproperties of media that are not accessible to traditional energysources used for transmission imaging techniques. For example, oneflourishing application of imaging in scattering media is in the fieldof optical tomography. Optical tomography permits the use of nearinfrared energy as an imaging source. Near infrared energy is highlyscattered by human tissue and is therefore an unsuitable source fortransmission imaging in human tissue. However, these properties make ita superior imaging source for scattering imaging techniques. The abilityto use near infrared energy as an imaging source is of particularinterest in clinical medicine because it is exceptionally responsive toblood volume and blood oxygenation levels, thus having great potentialfor detecting cardiovascular disease, tumors and other disease states.

A common approach for the reconstruction of an image of thecross-sectional properties of a scattering medium is to solve aperturbation equation based on the radiation transport equation. Theradiation transport equation is a mathematical expression describing thepropagation of energy through a scattering medium. The perturbationformulation relates the difference between coefficient values of thetrue target and a specified reference medium, weighted by aproportionality coefficient whose value depends on, among other things,the source/detector configuration and the optical properties of themedium. In practice, tomographic measurements consider some array ofmeasurement data, thus forming a system of linear equations having theformu−u _(r) =δu=W _(r) δx,  (1)where δu is the vector of differences between a set of measured lightintensities (u) and those predicted for a selected reference medium(u_(r)), W_(r) is the Jacobian operator, and δx is theposition-dependent difference between one or more optical properties ofthe target and reference media (i.e., a change in absorption coefficientδμ_(a), a change in the reduced scattering coefficient, μ′_(s), or, inthe diffusion approximation, the diffusion coefficient δD, whereD=1/[3(μ_(a)+μ′_(s))]). The operator, referred to as the weight matrix,has coefficient values that physically represent the fractional changein light intensity at the surface caused by an incremental change in theoptical properties at a specified point in the medium. Mathematicallythis is represented by the partial differential operator ∂u_(i)/∂x_(j),where i is related to the i^(th) source/detector pairs at the surface ofthe medium, and j to the j^(th) pixel or element in the medium.

Although the perturbation equation in Eq. (1) can be solved using any ofa number of available inversion schemes, practical experience has shownthat the accuracy and reliability of the results obtained can varygreatly due to uncertainties and errors associated with the quality ofthe measurement data, inaccuracies in the physical model describinglight propagation in tissue, specification of an insufficiently accuratereference state, the existence of an inherently underdetermined statecaused by insufficiently dense measurement sets, weak spatial gradientsin the weight function, and so forth.

In practice, a matter of considerable concern is the accuracy with whichthe reference medium can be chosen. An accurate reference is one thatclosely matches the external geometry of the target medium, has the samesize, nearly the same internal composition, and for which the locationsof the measuring probes and their efficiency coincide well with thoseused in the actual measurements. While such conditions may be easily metin numerical and perhaps laboratory phantom studies, they represent amuch greater challenge in the case of tissue studies. Confoundingfactors include the plasticity of tissue (it deforms upon probecontact), its mainly arbitrary external geometry and internalcomposition and the considerable uncertainty stemming from the expectedvariable coupling efficiency of light at the tissue surface. Theinfluence of these uncertainties can be appreciated when it isrecognized that the input data vector for the standard perturbationformulation (i.e., Eq. (1)) is actually the difference between ameasured and a computed quantity. This vector contains informationregarding the subsurface properties of the target medium that, inprinciple, can be extracted provided an accurate reference medium isavailable.

In practice, however, there are two significant concerns that arefrequently encountered in experimental studies and are not easilyresolvable especially in the case of tissue studies. One concern is theexpected variable coupling efficiency of light entering and exitingtissue. Nonuniformity in the tissue surface, the presence of hair orother blemishes, its variable deformation upon contact with opticalfibers, the expected variable reactivity of the vasculature in thevicinity of the measuring probe all serve to limit the ability toaccurately determine the in-coupling and out-coupling efficiencies ofthe penetrating energy. Consideration of this issue is critical asvariations in the coupling efficiency will be interpreted by thereconstruction methods as variations in properties of the target mediumand can introduce gross distortions in the recovered images. Inprinciple, the noted concern can be minimized by adopting absolutecalibration schemes, however, in practice the variability in tissuesurface qualities will limit reliability and stability of these efforts.

A second concern stems from the underlying physics of energy transportin highly scattering media. One effect of scattering is to greatlyincrease the pathlength of the propagating energy. Small changes in theestimated absorption or scattering properties of the medium can,depending on the distance separating the source and detector, greatlyinfluence the density of emerging energy. This consideration hasimportant implications regarding the required accuracy by which thereference medium must be specified. In the context of perturbationformulations, the reference medium serves to provide estimates of thepredicted energy density as well as to provide the needed weightfunctions that serve as the imaging operators. The difficulty is thatthe computed reference intensity values are extremely dependent on theoptical coefficient values of the reference medium. Significantly, thisdependence is a nonlinear function of the distance between source anddetector. It follows that a small change in the optical properties ofthe reference medium may influence the value of the computed intensitydifferences (δu) by a relative amount that may be significantlydifferent for each source/detector pair, thereby altering theinformation content of the data vectors. This can lead to the recoveryof grossly corrupted images. Whereas, in principle, such effects may beovercome by use of recursive solutions to the perturbation equation(i.e., Newton-type updates), in practice this can require extensivecomputational efforts, especially in the case of 3D solutions. Moreover,it is well appreciated that such efforts to improve on first ordersolutions to the perturbation equation (e.g., Born or Rytov solutions),can fail if the initial estimate chosen for the reference medium isinsufficiently accurate.

One alternative to devising absolute calibration schemes is to devisemethodologies whose solutions are intrinsically less sensitive, orbetter still, do not require such information, but nevertheless arecapable of providing accurate descriptions of certain features of highlyscattering media. While a range of empirical methodologies can bedevised, it is desirable that they be broadly extendable withoutrequiring undue physical approximations, since these are generallyincompatible with model-based methods.

An approach previously adopted is to directly relate relative detectorreadings, obtained from comparison of detector values derived from twodifferent target media (usually media with and without the includedobject), to the weight matrix computed based on a previously assignedreference medium. R. L. Barbour, H. Graber, R. Aronson, and J. Lubowsky,“Model for 3-D optical imaging of tissue,” Int. Geosci. and RemoteSensing Symp., (IGARSS), 2, 1395–1399 (1990). While capable of producinggood quality images of internal structure of a target medium, the methodproved to have limited utility as it did not produce solutions havingphysical units, thereby rendering specific interpretation difficult, aswell as limiting efforts to compute recursive solutions.

For the forgoing reasons, there is a need for image reconstructiontechniques based on the detection of scattered energy that (1) do notrequire absolute calibration of, and absolute measurements by, thedetectors and other elements of the apparatus, (2) make the standardperturbation equation less susceptible to variations between boundaryconditions and properties of the reference medium and the target medium,and (3) produce solutions having physical units.

SUMMARY

The present invention satisfies these needs by providing a method forgenerating an image of a scattering medium using normalized relativemeasured intensity values a perturbation formulation based on theradiation transport equation.

It is an object of the present invention to provide a method for imagingthe properties of a scattering target medium using a modifiedperturbation equation. The method comprises generating a first datavector of measured data from a target and a second vector of measureddata from a target, normalizing the first and second vectors of measureddata and solving a modified perturbation equation for the unknownoptical properties of a target medium. The first and second vectors ofmeasured data are measures of energy emerging from the target.

It is a further aspect of this invention to obtain the first and secondsets of measured data from the same target, wherein the first set ofmeasured data is a set of data measured at an instant in time and thesecond set of measured data is a time average mean of a plurality offirst sets of measured data.

It is yet a further aspect of this invention to obtain the first andsecond sets of measured data from two different targets.

In another aspect of this invention the modified perturbation equationis a modified Rytov approximation. In another aspect of this inventionthe modified perturbation equation is a modified Born approximation.

BRIEF DESCRIPTION OF THE FIGURES

For a better understanding of the invention, together with the variousfeatures and advantages thereof, reference should be made to thefollowing detailed description of the preferred embodiment and to theaccompanying drawings, wherein:

FIG. 1 is a schematic illustration of an exemplary imaging system;

FIG. 2A is a cross-sectional image of the absorption coefficients of thetarget;

FIG. 2B is a cross-sectional image of the diffusion coefficients of thetarget;

FIG. 3 is a table illustrating a summary of the test cases explored;

FIG. 4A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 1 of FIG. 3;

FIG. 4B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 1 of FIG. 3;

FIG. 5A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 2 of FIG. 3;

FIG. 5B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 2 of FIG. 3;

FIG. 6A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 3 of FIG. 3;

FIG. 6B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 3 of FIG. 3;

FIG. 7A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 5 of FIG. 3;

FIG. 7B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 5 of FIG. 3;

FIG. 8A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 6 of FIG. 3;

FIG. 8B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 6 of FIG. 3;

FIG. 9A is a series of reconstructed cross-sectional absorption profileimages of the target obtained in test case 7 of FIG. 3;

FIG. 9B is a series of reconstructed cross-sectional diffusion profileimages of the target obtained in test case 7 of FIG. 3;

FIG. 10 is a table listing constant calibration errors corresponding toeach image of FIGS. 9A and 9B;

FIG. 11 is a table of data corresponding to the error, resolution andcontrast for the reconstructed images shown in FIGS. 4A and B;

FIG. 12 is a table of data corresponding to the error, resolution andcontrast for the reconstructed images shown in FIGS. 5A and B;

FIG. 13A is a series of reconstructed cross-sectional absorption profileimages of the target shown in FIG. 15 with varying concentrations ofhemoglobin;

FIG. 13B is a series of reconstructed cross-sectional diffusion profileimages of the target shown in FIG. 15 with varying concentrations ofhemoglobin;

FIG. 14A is a series of reconstructed cross-sectional absorption profileimages of the target shown in FIG. 15 with varying reference mediumproperties;

FIG. 14B is a series of reconstructed cross-sectional diffusion profileimages of the target shown in FIG. 15 with varying reference mediumproperties;

FIG. 15 is a schematic illustration of a phantom study;

FIG. 16A is a graph plotting the amplitude of normalized intensitiesused in the modified perturbation formulation corresponding to thereconstructed cross-sectional images shown in row 3 of FIG. 4;

FIG. 16B is a graph plotting the amplitude of normalized intensitiesused in the modified perturbation formulation corresponding to thereconstructed cross-sectional images shown in column 3 of FIG. 4;

FIG. 17A is a graph plotting the amplitude of the relative intensitiesused in the standard perturbation formulation corresponding to thereconstructed cross-sectional images shown in row 3 of FIG. 5;

FIG. 17B is a graph plotting the amplitude of the relative intensitiesused in the standard perturbation formulation corresponding to thereconstructed cross-sectional images shown in column 3 of FIG. 5;

FIG. 18 is a graph plotting the amplitude of the frequency spectrum ofthe Fourier transforms for data vectors computed using equation (1) andequation (3); and

FIG. 19 is a table illustrating the ratio of average contrasts ofreconstructed absorption and diffusion coefficients shown in FIGS. 4Aand B respectively.

DETAILED DESCRIPTION OF THE INVENTION

Imaging in a scattering medium relates to the methods and techniques ofimage the internal properties of a medium based on the detection ofscattered energy. The typical process for imaging of a scattering mediumcomprises: (1) selecting a reference medium having known boundaryconditions and optical properties which are substantially similar tothose of the intended target; (2) determining a weight matrix and anintensity of emerging energy exiting the reference medium at each of aplurality of source points for each of a plurality of detectors locatedaround the reference medium boundary, the determination being made byeither actual measurements or solution of the radiation transportequation; (3) measuring actual emerging energy intensities forcorresponding source and detector points on a target medium; and (4)solving the perturbation equation for the optical properties of thetarget based on the measured intensities of energy emerging from thetarget.

The present invention describes an improved methodology for imaging of ascattering medium using a modified form of the standard perturbationequation. The inventive modification of the standard perturbationequation is capable of (1) reducing the sensitivity of the perturbationequation to differences between the reference medium and target medium,(2) producing solutions to the perturbation equation having physicalunits, and (3) reducing the effect of variable detector efficiencieswithout the need for absolute calibration, while at the same timepreserving the ability to compute recursive solutions. Compared to thestandard perturbation approach, the described invention providesremarkable improvement in the quality of image reconstruction.

While the method of the present invention is applicable to known staticimaging techniques it is instrumental in the realization of practicaldynamic imaging of highly scattering media. There are three principalelements to practical dynamic imaging. The first element is thedevelopment of a fast, parallel, multi-channel acquisition system thatemploys geometrically adaptive measuring heads. This system is describedbriefly below and in further detail in the copending Barbour 4147PC1application. The second element is to evaluate the acquired tomographicdata using the modified perturbation methods of the present invention.The third element is to collect a time series of data and subject eitherthe time series of data or a time series of reconstructed images fromthe data to analysis using various linear and nonlinear time-seriesanalysis methods to extract dynamic information and isolated dynamicinformation. These methods are described in detail in the copendingBarbour 4147PC2 application.

Some of the methods, systems and experimental results described belowfocus on optical tomography of human tissue using wavelengths in thenear infrared region for the imaging source. However, as disclosedgenerally herein, it will be appreciated to those skilled in the artthat the invention is applicable to the use of essentially any energysource (e.g., electromagnetic, acoustic, and the like) on any scatteringmedium (e.g., body tissues, oceans, foggy atmospheres, earth strata,industrial materials) so long as diffusive type mechanisms are theprincipal means for energy transport through the medium.

System

Numerous imaging systems such as those disclosed in the previouslymentioned the Barbour '355 patent, the Barbour '322 patent and theBarbour 4147PC1 application have been developed for use in imaging of ascattering medium. A schematic illustration of an exemplary system isshown in FIG. 1. This system includes a computer 102, sources 104, 106,a source demultiplexer 108, an imaging head 110, detectors 112 and adata acquisition board 114.

A target 116 placed in the imaging head 110 is exposed to optical energyfrom sources 104, 106. The optical energy originating from sources 104,106, is combined by beam splitter 118 and is delivered to sourcedemultiplexer 108. The source demultiplexer 108 is controlled bycomputer 102 to direct the optical energy to source fibers 120sequentially.

Each source fiber 120 carries the optical energy from the demultiplexer108 to the imaging head 110, where the optical energy is directed intothe target 116. The imaging head 110 contains a plurality of sourcefibers 120 and detector fibers 122 for transmitting and receiving lightenergy, respectively. Each source fiber 120 forms a source-detector pairwith each detector fiber 122 in the imaging head 110 to create aplurality of source-detector pairs. The optical energy entering thetarget 116 at one location is scattered and may emerge at any locationaround the target 116. The emerging optical energy is collected bydetector fibers 122 mounted in the imaging head 110.

The detector fibers 122 carry the emerging energy to detectors 112, suchas photodiodes or a CCD array, that measure the intensity of the opticalenergy and deliver a corresponding signal to a data acquisition board114. The data acquisition board 114, in turn, delivers the data tocomputer 102.

This imaging process is repeated so as to deliver optical energy to eachof the source fibers sequentially, a measurement being obtained fordetected emerging energy at each detector for each emitting sourcefiber. This process may continue over a period of time with the computer102 storing the data for reconstruction of one or more images.Additionally, the system may include two or more imaging heads forcomparing one target to another. The computer 102 reconstructs an imagerepresentative of the internal optical properties of the target bysolving a perturbation equation. It will be appreciated by those skilledin the art that more than one computer can be used to increase datahandling and image processing speeds.

The Standard Perturbation Formulation

As discussed above, reconstruction of a cross section image of theabsorption and/or scattering properties of the target medium is based onthe solution of a perturbation formulation of the radiation transportequation. The perturbation method assumes that the composition of theunknown target medium deviates only by a small amount from a knownreference medium. This reduces a highly non-linear problem to one thatis linear with respect to the difference in absorption and scatteringproperties between the target medium under investigation and thereference medium. The resulting standard perturbation equation has thefollowing forms:u−u _(r) =δu=W _(r) δx,  (1)In equation (1), δu is a vector of source-detector pair intensitydifferences between the measured target medium and the known referencemedium (i.e., δu=u−u_(r)). W is the weight matrix describing theinfluence that each volume element (“voxel”) of the reference medium hason energy traveling from each source to each detector, for allsource-detector pairs. The volume elements are formed by dividing aslice of the reference medium into an imaginary grid of contiguous,non-overlapping pieces. Physically, the weight matrix contains the firstorder partial derivatives of the detector responses with respect to theoptical coefficients of each volume element of the reference medium. δxis the vector of differences between the known optical properties (e.g.,absorption and scattering coefficients) of each volume element of thereference medium and the corresponding unknown optical properties ofeach volume element of the target medium.

This standard perturbation equation assumes (1) use of absolute detectormeasurements for u, and (2) that the boundary conditions and opticalproperties of the reference do not vary significantly from the target.Both of these factors are problematic in practice.

The Modified Perturbation Formulation

The present invention modifies the standard perturbation equation byreplacing δu with a proportionate relative difference between twomeasured values multiplied by a reference term of the required units asset forth in the equation (2) below:

$\begin{matrix}{\left( {\delta\; I_{r}} \right)_{i} = {\left\lbrack \frac{I_{i} - \left( I_{0} \right)_{i}}{\left( I_{0} \right)_{i}} \right\rbrack\left( I_{r} \right)_{i}}} & (2)\end{matrix}$where i is the source/detector pair index. In equation (2), I_(r) is thecomputed detector reading corresponding to a source-detector pair of aselected reference medium, and I and I₀ represent two data measurementsfor a corresponding source-detector pair on one or more targets (e.g.,background vs. target, or time-averaged mean vs. a specific time point,etc.). The resultant term δI_(r) therefore represents a relativedifference between two sets of measured data that is then mapped onto areference medium. Careful examination reveals that this modification hasimportant attributes that limit the effects of modeling errors andminimize ill-conditioning of the inverse problem while retaining thecorrect units in the solution.

The corresponding perturbation equation using this modified term is:W _(r) ·δx=δI _(r)  (3)In equation (3) W_(r) and δx are the same as W_(r) and δu in equation(1). Referring to equations (2) and (3), it can be seen that in thelimit, when I_(r) equals to I₀, this equation reduces to the standardperturbation formulation shown in equation (1). This formulationrepresents the Born approximation formulation of the modifiedperturbation equation. A similar expression may be written for the Rytovapproximation in the following form:

$\begin{matrix}{{{\left( {\delta\; I^{\prime}} \right)_{i} = {\ln\frac{I_{i}}{\left( I_{0} \right)_{i}}}};}{{\left( W_{r}^{\prime} \right)_{ij} = \frac{\left( W_{r} \right)_{ij}}{\left( I_{r} \right)_{i}}};}{{\delta\; I^{\prime}} = {W_{r}^{\prime}\delta\; x}}} & (4)\end{matrix}$The inventive operation accomplished by equation (2) is to preserve theinformation content of a measured proportionate relative data vectorsobtained from the target medium and to yield data vectors having thecorrect physical units. Apart from simplifying measurement requirements,the method represented by equations (3) and (4) also reducessusceptibility of the perturbation equation to boundary and opticalproperty variation between the target and the reference medium.Additionally, iterative solutions of equations (3) and (4) can be easilyimplemented. In principle, the techniques used in the modifiedperturbation equation, referred to as the normalized difference method(NDM), can be used to evaluate any measured proportionate relativequantity.Experimental Validation

The following discussion presents results validating the methods andadvantages of the present inventions. These examples are presentedmerely as an illustration of the benefits of applying the NDM approachfor the analysis of relative measures from highly scattering media.

Forward Model and Data Acquisition Geometry

For any intended target the perturbation equation is generated for areference medium having boundary conditions and optical propertiessubstantially similar to the target. The perturbation equation modelsthe energy propagation, e.g. light, in the reference medium as adiffusion process. For a domain Ω having a boundary ∂Ω, this isrepresented by the expression:∇·[D(r)∇u(r)]−μ_(a)(r)u(r)=δ(r−r _(s)),rεΩ  (5)where u(r) is the photon density at position r, r_(s) is the position ofa DC point source, and D(r) and μ_(a)(r) are the position-dependentdiffusion and absorption coefficients, respectively. The diffusioncoefficient is defined as:D=1/[3(μ_(a) [r]+μ _(s) [r])]  (6)where μ′_(s)[r] is the reduced scattering coefficient. The photondensity values at the detectors, i.e., the calculated energy intensityemerging from the reference medium at each detector, were computed byapplying Dirichlet boundary conditions on an extrapolated boundary.Depending on the target medium to be explored, sources and detectors forthe reference are positioned 1 to 2 transport mean free pathlengthswithin the boundary of the reference medium.

Solutions to the diffusion equation may be computed by any known means,such as by the KASKADE adaptive finite element method. R. Beck, R.Erdmann and R. Roitzsch, “Kaskade 3.0—An object-oriented adaptive finiteelement code,” Technical report TR 95-4, Konrad-Zuse-Zentrum furInformationstechnik, Berlin (1995). This is a publicly available codesuitable for the solution of partial differential equations in one, twoor three dimensions using adaptive finite element techniques. The codecan be readily modified to permit solutions to the diffusion equationusing a point source. Mesh generation may be by any known method, suchas the Delaunay tessellation algorithm originally proposed by Watson. D.F. Watson, “Computing the n-dimensional Delaunay tessellation withapplications to Voronoi polytopes”, Computer Journal, 24, 167–172(1981).

The perturbation equation is specific to the boundary conditions andoptical properties of the reference medium, including the orientation ofthe source-detector pairs in relation to one another and the referencemedium. These conditions and properties are preferably nearly identicalto the target. For example, in the experiments discussed below, theperturbation equation was generated based on an imaging system havingsix sources and eighteen detectors per source (108 source-detectorpairs) with the sources equally spaced at 60 degree intervals around theboundary of the medium and the detectors equally spaced at 20 degreeintervals.

Inverse Algorithm

As described above, in the present invention relative intensities aremeasured for all source-detector pairs using any known imaging system.Image recovery is then achieved using known methods, such as conjugategradient descent (CGD), or simultaneous algebraic reconstructiontechniques (SART), to solve the modified perturbation equation for theabsorption and scattering properties of the target. J. Chang, H. L.Graber, R. L. Barbour and R. Aronson, “Recovery of optical cross-sectionperturbations in dense-scattering media by transport-theory-basedimaging operators and steady-state simulate data”, Appl. Opt. 35,3963–3978, (1996) (the disclosure of which is incorporated herein byreference). For example, the experimental results discussed below wereachieved using a CGD solver with and without matrix resealing. Inaddition, a weight matrix rescaling (WMR) technique may be used toimprove the ill-conditioning of the weight matrix. The effect ofresealing the weight matrix is to make it more uniform. Two rescalingcriteria can be applied for this purpose: (1) rescaling the maximum ofeach column to 1; or (2) rescaling the average of each column to 1. Inthe experimental results below, when WMR was used, criterion 1 wasapplied for image recovery. The solution to the modified perturbationequation provides a relative measure of the difference between thecross-sectional optical properties of a target during the first andsecond measurements I and I₀. The values from this solution are used togenerate cross-sectional images representative of the target's internaloptical properties.

Test Cases Explored

The following discussion presents results obtained for seven test casescomparing image reconstruction using the known standard perturbationformulation with the modified perturbation formulation of the presentinvention. These examples are presented merely as an illustration of thebenefits of the modified perturbation method of the present invention.

The reconstruction results presented in the test case are limited tosolution of the first order Born approximation. The coefficient valuesfor absorption and diffusion coefficients were computed simultaneously.For each case tested, measures of error, contrast accuracy andresolution were also computed. These are defined as follows:

-   -   i. Image error is the relative root mean square error (RMSE)

$\begin{matrix}{{RMSE} = \sqrt{\frac{\sum\limits_{j = 1}^{M}\;\left( {x_{j} - a_{j}} \right)^{2}}{\sum\limits_{j = 1}^{M}\;\left( a_{j} \right)^{2}}}} & (7)\end{matrix}$where a_(j) and x_(j) are the actual and reconstructed values of theoptical coefficient, and M is the number of volume elements used forreconstruction.

-   -   ii. Image contrast accuracy was determined by computing the mean        value of the recovered perturbation coefficient along the        transect bisecting the two objects.    -   iii. Resolution was measured by computing the mean value of the        full-width half-maximum of the two reconstructed objects along        the transect bisecting the inclusions.

FIGS. 2A and 2B show the cross-sectional geometry and absorption (FIG.2A) and diffusion (FIG. 2B) coefficient profiles of the target mediumexplored. The target medium is 8 cm in diameter and has two includedobjects each 1 cm in diameter and separated by 3 cm symmetrically aboutthe center of a homogeneous background medium. Optical properties of thebackground and included objects are 0.04 and 0.02 cm⁻¹ for μ_(a)(absorption coefficient), and 10 and 5 cm⁻¹ for μ_(s) (scatteringcoefficient), respectively.

The table of FIG. 3 lists the various test cases explored. The symbol“V” indicates that the parameter was varied, “C” indicates that theparameter was held constant, and “/” indicates that the parameter wasnot considered. The test cases allowed at least partial isolation of theeffect of variations in each of the input parameters on the resultantimage for different reconstruction schemes and perturbationformulations. This testing permitted exploration of the dependence ofill-conditioning on the different input parameters that influence theaccuracy and stability of the image reconstruction.

Test cases 1 and 2 examined the general case where the reference mediumis based only on an estimate of the background optical properties of thetarget medium. The estimated properties were varied over a broad range,ranging from 0.0 cm⁻¹ to 0.3 cm⁻¹ in μ_(a) and from 3 cm⁻¹ to 30 cm⁻¹ inμ_(s). For purposes of comparison, test cases 3 through 7 explored thedependence of image quality on the varied parameters using the standardperturbation formulation.

Test cases 3 and 4 mainly mirror conditions explored in cases 1 and 2with the exceptions that the standard perturbation formulation wasevaluated, and a narrower range of coefficient values was considered forthe reference medium. Here the general case is also considered whereonly an estimate of the background optical properties of the targetmedium is available. The range of values for the optical propertiesexplored were from 0.02 cm⁻¹ to 0.08 cm⁻¹ in μ_(a), and from 5 cm⁻¹ to15 cm⁻¹ in μ_(s).

Cases 5 and 6 consider the special situation where prior knowledge ofthe background properties of the target medium is known. The parametersvaried were W_(r) and I_(r), referred to as W_(b) and I_(b),respectively. The range of optical properties varied for test 5 is sameas in case 1. For case 6, the range of optical properties varied is thesame as in case 3. Test case 7 explores the effect of a constantcalibration error in measurement, and assumes prior knowledge of thebackground properties of the target medium.

Results

Data shown in FIGS. 4 through 9 illustrates the influence that thevaried parameters listed in FIG. 3 have on the reconstruction resultsderived from a first-order Born approximation using the standard andmodified perturbation formulations. The results presented are listed ina matrix format. The value of the absorption and scattering coefficientfor the reference medium is fixed for each row and column, respectively.Varied is the value of these parameters along the orthogonal direction.Shown in the figures are the reconstruction profiles for all test casesexplored except case 4, whose findings are reported in the text.

Qualitative Analysis

Data in FIGS. 4 and 5 show the quality of reconstructed images obtainedusing equation (3). FIGS. 4A and 5A illustrate the computed absorptionmaps, while FIGS. 4B and 5B show the computed diffusion maps. Comparedto the original profiles shown in FIG. 2, these findings illustrate thatqualitatively accurate results, revealing the two-object structure, areobtained over a broad range of values for the selected reference medium.Artifact dominated results are limited to cases in the lower rightcorner of the matrix. These correspond to those reference media havingabsorption and scattering coefficient values significantly greater thanthe background of the target medium. Quantitative analysis of theseresults is presented in the next section. Comparison of images shown inFIGS. 4 and 5 reveals that the matrix rescaling method is capable ofproviding a higher resolution image, though over a reduced range ofvalues for the reference medium. For example, comparison of resultsreveals that the matrix resealing method yields only artifacts in theabsorption map for non-absorbing reference media, while under the sameconditions the diffusion coefficient maps reveal two completely resolvedobjects. In the absence of matrix rescaling, both coefficient mapsreveal the presence of the included objects, though with reduced edgeresolution and more artifacts in the diffusion map. The addedimprovement using matrix rescaling, however, is achieved under thelimiting conditions of a range constraint (positive for D and negativefor μ_(a)). Overall, the range of values for the reference medium'soptical properties for which qualitatively accurate maps are obtainedfor both inverse methods is much greater than previously reported. S. R.Arridge, M. Schweiger, “Sensitivity to prior knowledge in opticaltomographic reconstruction”, in Proc. Optical tomography, photonmigration and spectroscopy of tissue and model media: Theory, humanstudies, and instrumentation, SPIE, 2389, 378–388, (1995) (thedisclosure of which is incorporated herein by reference).

Shown in FIG. 6 are results from case 3 evaluated using equation (1).Compared to results shown in FIGS. 4 and 5, a more limited range ofvalues of optical properties for the reference medium were examined,since outside of this range, only artifact was recovered. Even withinthe explored range, significant instability was observed for relativelysmall variations in the reference medium. This sensitivity indicates astate of ill-conditioning that is alleviated using the modifiedperturbation formulation (equation 3). Not shown are results from case 4using the matrix rescaling method. In case 4, even greater instabilitywas observed than in the case using CGD only.

The parameters varied in the above figures include both the computedreference intensity and the weight matrix. This is the general casewhere both quantities can only be estimated and are computed from aspecified reference medium.

Results shown in FIGS. 7 through 9 explore the special cases whereerrors occur only in one parameter. Data in FIG. 7 illustrates theinfluence of errors in the estimated weight matrix. Assumed is priorknowledge of the reference detector intensity values for the backgroundmedium. In practice this would correspond to situations where ameasurement was made in the presence and absence of an included object.Inspection reveals that qualitatively accurate results are obtained overa significantly broader range of reference values than those presentedin FIG. 7. This finding suggests that the principal origin of theill-conditioning of the inverse problem is associated with errors in theestimated reference intensity.

Results shown in FIG. 8 explore this possibility directly. In thissituation, we assume the unlikely case where accurate prior knowledge ofthe weight matrix for the target medium is available. Comparison of theresults in FIG. 8 to the results in FIG. 6 illustrates some improvementin the range of reference media yielding qualitatively accurate results.However, this range is small compared to the situations tested in FIG. 8for the standard perturbation formulation and for the modifiedperturbation formulation in FIGS. 4 and 5.

Finally, in FIG. 9, we test the case where accurate prior knowledge forthe weight matrix and reference detector values are available and aconstant error of measurement is introduced. The error added to themeasured detector reading, I, varies from −50 to 900%. Variations inconstant calibration error corresponding to FIG. 8 are shown in FIG. 9.Under these conditions, the results show that a constant calibrationerror does not significantly affect the qualitative accuracy of thecomputed images. It is worth noting that errors of this type will notexist using the modified perturbation formulation.

Quantitative Analysis

A quantitative analysis of the image data was made by computing measuresof error, resolution and contrast. Results shown in FIGS. 11 and 12 arethe corresponding values computed from data shown in FIGS. 4 and 5. Theformat of the data is the same as in FIG. 4 (i.e., data in the rows arederived for a fixed value of absorption, while column data is derivedfor a fixed value of scattering).

Inspection of FIGS. 11 and 12 reveals, not surprisingly, that the lowestimage RMSE is achieved when the selected reference medium matches thebackground optical properties of the target medium. Consistent with thisfinding is that nearly equivalent error values were obtained for thosemaps in which one of the fixed coefficient values for the referencemedium matched that of the background. Interestingly, in the absence ofWMR, whereas these conditions produced the lowest total error values forthe image map, improved accuracy of object contrast was obtained usingreference media generally having reduced scattering and increasedabsorption values compared to those of the background. An exception tothis was the case where the background absorption level in the referencemedium was reduced while scattering value matched the object. With WMR,the best accuracy for object coefficient values was achieved when thereference medium matched the coefficient values of the backgroundmedium. This is not unexpected given the imposed constraint. Overallthere is evidence of a trade-off between artifact levels and accuracy inobject contrast. Potentially significant is the observation of severalinstances where considerably enhanced object contrast is seen withoutundue degradation of image quality. This is observed with either inversealgorithm, though the trends are somewhat different. In the absence ofWMR, enhanced contrast for both absorption and diffusion maps are seenusing reference media whose absorption and scattering coefficient valuesare lower than the background. Enhanced contrast is also seen in thecase of a non-absorbing reference medium, although increased artifactsare present. With WMR, improved object contrast is also seen withreference media having reduced scattering, but the trend in contrastenhancement is opposite for the different coefficients and dependsstrongly on the value of the absorption coefficient. Increasing theabsorption coefficient value for the reference medium increases theobject contrast value for absorption while reducing the contrast for thediffusion coefficient. Inspection of results reported for edgeresolution reveal an under and over estimate of object diameter forimages computed with and without the WMR, respectively. Whereas it isoften best to avoid errors in edge resolution, the resolution of theimages obtained using the WMR method is striking.

Experimental Validation of Modified Perturbation (NDM) Formulation

Experimental verification demonstrating that the modified perturbationformulation is capable of resolving internal structure of a densescattering medium is given in FIGS. 13 and 14. Tomographic measurementswere performed at 780 nm using the IRIS imaging system previouslydescribed. R. L. Barbour, R. Andronica, Q. Sha, H. L. Graber, and I.Soller, “Development and evaluation of the IRIS-OPTIsoanner, ageneral-purpose optical tomographic imaging system.” In Advance inOptical Imaging and Photon Migration, J. G. Fujimoto et al, ed., Vol. 21of OSA Proceeding Series, pp. 251–255, (1998) (the disclosures of whichare incorporated herein by reference). The target medium was a latexlaboratory glove filled with 2% (v/v) Intralipid suspended from a holderin a pendant position. Added to the glove were two 1 cm diameter plastictubes filled with varying concentrations of hemoglobin (Hb) in theamount of 5 μm, 10 μm, 20 μm, and 40 μm. A cross section of the phantomset-up is shown in FIG. 15. The pass-through diameter of the IRISimaging head was closed until gentle contact with the glove wasachieved. The diameter of the glove in the measurement plane was 6.7 cm.Above and below this plane, the glove assumed an arbitrary geometry.Tomographic measurements were performed using the same measurementgeometry described for the numerical studies. Optical measurements wereperformed in the presence and absence of the included objects from whichthe relative intensity values were derived. The resultant data vectorswere then evaluated by equation (3) using a regularized CGD methodwithout weight matrix resealing. The coefficient values for thereference medium used were varied from 0.01 to 0.04 cm⁻¹ in μ_(a) and 10to 20 cm⁻¹ in μ_(s).

Reconstructed μ_(a) and D maps for each experiment with differentconcentrations of hemoglobin in the two tubes using a specific referencemedium are shown in FIG. 13. Inspection reveals two well-resolvedobjects whose contrast in both coefficients increases with increasingconcentration of added absorber. Quantitatively, the dependence of imagecontrast on absorber concentration is less than expected (δμ_(a)recovered=0.005 vs. δμ_(a) actual=0.015 cm⁻¹) indicating perhaps eithera self-shielding effect, limits of a first order Born solution or both.Results shown in FIG. 14 demonstrate that images of similar quality arederivable over a range of reference media, a result consistent with theabove described numerical studies.

DISCUSSION AND CONCLUSIONS

The present invention describes and evaluates a new formulation for theinverse problem for imaging in highly scattering media. Motivating thisdevelopment has been an appreciation of the expected limits imposed bypractical measurements, especially as it relates to the dependence ofimage accuracy on instrument calibration and the ability to specify anaccurate reference medium. This concern arises because many of theanticipated clinical applications will require some level of accuracy inthe computed coefficient values. An accurate solution will require anexplicit accounting of various factors intrinsic to the detector (e.g.,quantum efficiency, acceptance angle etc.), as well as features specificto the target. This includes, in particular, the efficiency of contactwith the target medium by the detector or intervening optical fibersthat deliver and collect the optical signal. This is necessary becauseall model-based imaging schemes proposed for imaging in highlyscattering media assume equivalency in detector efficiency for measuredand computed values. Failure to account for such variables willintroduce error whose magnitude and stability could vary considerablydepending on the specifics of the measuring device and target medium.

In principle, these uncertainties can be taken into account, althoughnot without considerable effort for instrument design and addedcomplexity of the forward modeling code. Generally speaking, suchlimitations are widely appreciated by many, both in this and otherimaging communities. The goal in identifying practical schemes is todevise strategies that are mainly insensitive to such uncertainties.Often a desirable starting point is to employ schemes that provideuseful information based on some type of relative measurement.Previously, we described a back projection scheme that evaluatedrelative detector data. R. L. Barbour, H. Graber, R. Aronson, and J.Lubowski, “Model for 3-D optical imaging of tissue,” Int. Geosci andRemote Sensing Symp., (IGARSS), 2, 1395–1399 (1990) (the disclosures ofwhich are incorporated herein by reference). H. L. Graber, J. Chang, J.Lubowsky, R. Aronson and R. L. Barbour, “Near infrared absorptionimaging in dense scattering media by steady-state diffusion tomography”,in Proc. Photon migration and imaging in random media and tissues”,1888, 372–386, (1993) (the disclosures of which are incorporated hereinby reference). This formulation employed model-based imaging operators,but produced solutions lacking physical units. The lack of physicalunits (1) makes specific interpretation difficult, especially shouldmulti-wavelength measurements be considered, and (2) makes efforts tocompute iterative updates difficult. However, this scheme however showedthat in all cases tested, high contrast images having excellent edgedetection and object location were achievable. R. L. Barbour, H. Graber,R. Aronson, and J. Lubowski, “Model for 3-D optical imaging of tissue,”Int. Geosci and Remote Sensing Symp., (IGARSS), 2, 1395–1399 (1990) (thedisclosures of which are incorporated herein by reference). H. L.Graber, J. Chang, J. Lubowsky, R. Aronson and R. L. Barbour, “Nearinfrared absorption imaging in dense scattering media by steady-statediffusion tomography”, in Proc. Photon migration and imaging in randommedia and tissues “, 1888, 372–386, (1993) (the disclosures of which areincorporated herein by reference). H. L. Graber, J. Chang and R. L.Barbour “Imaging of multiple targets in dense scattering media”, inProc. Experimental and numerical methods for solving ill-posed inverseproblems: Medical and non-medical applications, SPIE, 2570, 219–234,(1995) (the disclosures of which are incorporated herein by reference).

In subsequent studies we have identified that the expression evaluatedusing this method has a functional form of:

$\begin{matrix}{{\left( \frac{\delta\; I}{I_{0}} \right)_{i}\left( {\sum\limits_{j = 1}^{N}\;\left( W_{r} \right)_{ij}} \right)} = {\sum\limits_{j = 1}^{N}\;\left( {\left( W_{r} \right)_{ij}\left( {\delta\; x} \right)_{j}} \right)}} & (8)\end{matrix}$

-   -   for each source-detector pair where, i is the source-detector        pair number, j is the element number and N is the number of        elements. It is apparent that this expression is similar to the        following equation, which is equation (3) with a different form,

$\begin{matrix}{{\left( \frac{\delta\; I}{I_{0}} \right)_{i}\left( I_{r} \right)_{i}} = {\sum\limits_{j = 1}^{N}\;\left( {\left( W_{r} \right)_{ij}\left( {\delta\; x} \right)_{j}} \right)}} & (9)\end{matrix}$

While it is evident the two expressions are not equivalent, a morecareful examination reveals that the different quantities on the lefthand side of Eqs. (8) and (9) (i.e., the sum of weights and thereference intensity for the i^(th) source-detector pair) are closelyrelated. The influence that different forms of the data vector have onimage recovery is discussed subsequently.

Solution of the perturbation formulation requires specification of threeinput data sets. Two quantities, I_(r) and W_(r), are typically computedfrom a specified reference medium, and the third quantity is themeasured response, I. Results presented in FIGS. 6 through 9 haveexplored the influence that errors in each of these quantities have onthe quality of computed reconstructions for the original perturbationformulation (equation (1)). Most sensitive was the case where bothquantities I_(r) and W_(r) are in error (c.f., test case 3 and FIG. 6).This suggests that the origin of excessive ill-conditioning can betraced to one or both of the quantities I_(r) and W_(r). When priorknowledge W_(r) is assumed, errors in I_(r) (test case 6) still producedhighly unstable solutions (FIG. 8). On the other hand, when priorknowledge of the correct values of I_(r) is assumed (test case 5), thesensitivity to errors in W_(r) is much less (FIG. 7), at leastqualitatively. This strongly suggests that the principal cause ofinstability can be traced to errors in I_(r). It also indicates thatshould error free measurement data be available, significant instabilityin the solution domain would persist. The influence of systematic errorin measurement was explored in test case 7. Results in FIG. 9 showedthat under conditions where prior knowledge of I_(r) and W_(r) for thebackground is available, qualitatively accurate solutions could beobtained even in the presence of 900% error. In the absence of thisprior knowledge, additional error due to this quantity further degradedimage quality even in those cases where the selected reference mediumdiffered only minimally from the background.

While we have demonstrated that the described reconstruction procedurecan significantly stabilize the computed reconstructions to errors inthe reference medium, and that the principal origin of solutioninstability in its absence is the result of to errors in I_(r), it isuseful to gain further insight as to why this should occur. An importantdifference between the two perturbation formulations is that equation(1) computes the difference between two exponentially attenuatedquantities, while equation (3) performs a linear operation on anexponentially attenuated quantity. Because of the non-linearrelationship between the medium coefficient values and surface detectorresponses, small errors in the former (the selected reference medium)can lead to large errors in the latter (the computed intensity or weightassociated with the selected reference medium). Moreover, because therelationship is non-linear, such errors can be expected to effectivelydistort the information content of the resultant data vector. Theoccurrence of such distortion is shown in FIGS. 16 and 17. FIG. 16 showsthe angle dependence of the relative detector response (i.e., the datavector) for a source bisecting the two inclusions computed usingreference media corresponding to row 3 and column 3 of FIG. 5.Inspection of these plots clearly reveal a bimodal attenuation profileindicating the presence of two buried objects, a finding consistent withits actual structure. In contrast, this structure is almost completelyabsent from results derived using the original perturbation formulation(see FIG. 17) even though the variation in range of the reference mediumis much less than that used in FIG. 16. This difference is most evidentin results given in FIG. 18 that shows the amplitude of the frequencyspectrum of the corresponding Fourier transforms. Comparison revealsthat in the case where the selected reference medium has an error onlyin μ_(a) of 0.02 cm⁻¹, this produces an error of approximately 5 ordersof magnitude in the amplitude of the frequency spectrum. This errorgrows to approximately eight orders of magnitude for the case where theselected reference medium has an error in μ_(s) of only 5 cm⁻¹. Giventhe magnitude of these errors, the extreme sensitivity of thereconstruction results obtained from the original perturbation equationis evident.

When we examined the computed data vector for the formulation based onequation (7), we observed a pattern similar to that shown in FIG. 16,but with much higher amplitude values. It is this close relationshipthat we believe accounts for why previous reconstruction results derivedusing equation (8) also provided stable and qualitatively accurate maps,even though the solutions lacked features important for specificinterpretation and iterative updates. R. L. Barbour, H. Graber, R.Aronson, and J. Lubowski, “Model for 3-D optical imaging of tissue,”Int. Geosci and Remote Sensing Symp., (IGARSS), 2, 1395–1399 (1990) (thedisclosures of which are incorporated herein by reference). H. L.Graber, J. Chang, J. Lubowsky, R. Aronson and R. L. Barbour, “Nearinfrared absorption imaging in dense scattering media by steady-statediffusion tomography”, in Proc. Photon migration and imaging in randommedia and tissues”, 1888, 372–386, (1993) (the disclosures of which areincorporated herein by reference). Y. L. Graber, J. Chang and R. L.Barbour “Imaging of multiple targets in dense scattering media”, inProc. Experimental and numerical methods for solving ill-posed inverseproblems: Medical and non-medical applications, SPIE, 2570, 219–234,(1995) (the disclosures of which are incorporated herein by reference).

A further advantage of the current method compared to the previouslydescribed SART-Type algorithm is that we are able to directly evaluateintensity difference values for which no mismatches exist between thecomputed intensity values, I_(r), and the Jacobian matrix, W_(r). Y. L.Graber, J. Chang and R. L. Barbour “Imaging of multiple targets in densescattering media”, in Proc. Experimental and numerical methods forsolving ill-posed inverse problems: Medical and non-medicalapplications, SPIE, 2570, 219–234, (1995) (the disclosures of which areincorporated herein by reference). Not only does this reduce systematicerrors, but produces solutions that can be updated by iterative methods.As reported here, error analysis studies have demonstrated that thecurrent methodology produces remarkable stable solutions havingexcellent qualitative accuracy.

An important goal shared by many investigators in the biomedical opticsfield is the capacity to accurately quantify variations in opticalcoefficients in tissue. One consequence of adopting the describedreconstruction procedure is that the derived solution will not convergeto the actual value even with error free measurement data and non-linearupdates. Instead, as is evident from results shown in FIGS. 4 and 5,they will converge to a solution that is proportional to the true valuethroughout the image map.

Examination of the recovered values indicates that the value of thisproportionality coefficient depends strongly on the selected referencemedium. Interestingly, though, as shown in FIG. 19, the ratio of thecomputed coefficient values is nearly constant over a broad rangeindicating that the relative error in the coefficients is itself aconstant. In fact, for the examples studied, we find that(δμ_(a)/δD)_(Reconstructed)≈(δμ_(a)/δD)^(1/2) _(Original) over most ofthe range of reference values explored. In other examples (results notshown), we have explored this relationship for a variety of perturbationvalues with a similar range of reference media. In all cases, a constanterror in the ratio of the derived coefficients was obtained. The valueof the proportionality constant, however, varied depending on themagnitude and direction of the perturbation, but remained within arelatively small range.

While it may be that case for some studies the restrictions imposed bythe described method may prove limiting, it is anticipated that thereare many practical situations where measurement of relative changes arenonetheless very useful. One example of special interest is in theevaluation of dynamic imaging data.

Finally, from a mathematical perspective, the described methodology hasthe effect of desensitizing the solution of boundary value problems tospecific features of the boundary. One practical consequence of this isthat, unlike the standard perturbation formulation, the currentformulation is less sensitive to a detailed knowledge of the externalboundary of tissue, a quantity not easily measured.

While the instant invention has considered using of DC sources, it isexpressly understood that the described methodology is extendable toother source conditions (e.g., time-harmonic and time-resolvedmeasurements) and other inverse formulations (e.g., iterative gradientdescent) and is extendable to other imaging modalities including,ultrasound, radioscintigraphic and impedance imaging.

Although illustrative embodiments have been described herein in detail,those skilled in the art will appreciate that variations may be madewithout departing from the spirit and scope of this invention. Moreover,unless otherwise specifically stated, the terms and expressions usedherein are terms of description and not terms of limitation, and are notintended to exclude any equivalents of the system and methods set forthin the following claims.

1. A computer-implemented method for imaging of the properties of atleast one scattering target medium, comprising: generating a firstvector of measured data and a second vector of measured data, the firstvector being indicative of energy emerging from the at least onescattering target medium, the second vector being indicative of energyemerging from the at least one scattering target medium, the emergingenergy substantially originating from at least one source directing theenergy into the at least one scattering target medium; normalizing adifference between the first and second vectors according to a ratio of:(a) a difference between the first and second vectors, to (b) the secondvector; solving a modified perturbation equation of a radiationtransport inverse problem for a relative change between a known propertyof a reference medium and a corresponding unknown property of the atleast one scattering target medium, wherein the modified perturbationequation relates the normalized difference and a vector of referencedata for the known reference medium to the relative change in theproperty, the vector of reference data being indicative of energyemerging from the known reference medium; and generating an image of theat least one scattering target medium, responsive to the solving.
 2. Thecomputer-implemented method of claim 1 wherein energy sources and energydetectors are configured in corresponding source-detector pairs for theat least one scattering target medium and the reference medium, and themodified perturbation equation has the following form:${\left( {\delta\; I_{r}} \right)_{i} = {\left\lbrack \frac{I_{i} - \left( I_{0} \right)_{i}}{\left( I_{0} \right)_{i}} \right\rbrack\left( I_{r} \right)_{i}}};\mspace{14mu}{and}$W_(r) ⋅ δ x = δ I_(r) where δx is a vector of the relative changesbetween the known property of the reference medium and the correspondingunknown property of the at least one scattering target medium, W_(r) isa weight matrix describing the influence that each of a plurality ofvolume elements of the reference medium has on energy emerging at apoint on the reference medium; I_(r), is the vector of reference data;(I_(r))_(i) is an element of I_(r) for an ith source-detector pair ofthe reference medium; I_(i) is an element of the first vector for an ithsource-detector pair of the at least one scattering target mediumcorresponding to the ith source-detector pair of the reference medium;I₀ is the second vector; (I₀)_(i) is an element of I₀ for the ithsource-detector pair of the at least one scattering target medium;(δI_(r))_(i) represents a change in the vector of reference data for theith source-detector pair of the reference medium; and δI_(r) representsa change in the vector of reference data for all source-detector pairsof the reference medium.
 3. The computer-implemented method of claim 1wherein the normalizing the difference between the first and secondvectors comprises determining a logarithm of ratio of the first vectorto the second vector.
 4. The computer-implemented method of claim 1wherein energy sources and energy detectors are configured incorresponding source-detector pairs for the at least one scatteringtarget medium and the reference medium, and the modified perturbationequation has the following form:${\left( {\delta\; I^{\prime}} \right)_{i} = {\ln\frac{I_{i}}{\left( I_{0} \right)_{i}}}};$${\left( W_{r}^{\prime} \right)_{ij} = \frac{\left( W_{r} \right)_{ij}}{\left( I_{r} \right)_{i}}};$δ I^(′) = W_(r)^(′)δ x where δx is a vector of the relative changesbetween the known property of the reference medium and the correspondingunknown property of the at least one scattering target medium, W′_(r) isa weight matrix describing the influence that each of a plurality ofvolume elements of the reference medium has on energy emerging at apoint on the reference medium; (W′_(r))_(ij) is an element of W′_(r) foran ith source-detector pair of the reference medium and a jth volumeelement of the reference medium; I_(r) is the vector of reference data;(I_(r))_(i) is an element of I_(r) for the ith source-detector pair ofthe reference medium; I_(i) is an element of the first vector for an ithsource-detector pair of the at least one scattering target mediumcorresponding to the ith source-detector pair of the reference medium;I₀ is the second vector; (I₀)_(i) is an element of I₀ for the ithsource-detector pair of the at least one scattering target medium; andδI′ is a vector representing a relative difference between the first andsecond vectors that is mapped onto the reference medium.
 5. Thecomputer-implemented method of claim 1 wherein the unknown propertywhich is solved for by the solving includes is at least one of anabsorption coefficient and a scattering coefficient.
 6. Thecomputer-implemented method of claim 1 wherein the first and secondvectors are obtained from one scattering target medium.
 7. Thecomputer-implemented method of claim 1 wherein the first and secondvectors are obtained from different scattering target mediums.
 8. Thecomputer-implemented method of claim 1 wherein the first and secondvectors are obtained at different time instants to provide dynamicimaging data of the at least one scattering target medium.
 9. Thecomputer-implemented method of claim 1 wherein the first vector isobtained at a first instant in time and the second vector represents atime averaged mean of a plurality of measurements.
 10. Thecomputer-implemented method of claim 1 further comprising generating across-sectional image representing the relative changes in the property.11. The computer-implemented method of claim 1, wherein solving themodified perturbation equation comprises mapping the normalizeddifference onto the reference medium.
 12. A system for imaging of theproperties of at least one scattering target medium, comprising: meansfor generating a first vector of measured data and a second vector ofmeasured data, the first vector being indicative of energy emerging fromthe at least one scattering target medium, the second vector beingindicative of energy emerging from at least one scattering targetmedium, the emerging energy substantially originating from at least onesource directing the energy into the at least one scattering targetmedium; means for normalizing a difference between the first and secondvectors according to a ratio of: (a) a difference between the first andsecond vectors, to (b) the second vector; means for solving a modifiedperturbation equation of a radiation transport inverse problem for arelative change between a known property of a reference medium and acorresponding unknown property of at least one scattering target medium,wherein the modified perturbation equation relates the normalizeddifference and a vector of reference data for the known reference mediumto the relative change in the property, the vector of reference databeing indicative of energy emerging from the known reference medium; andmeans for generating an image of the at least one scattering targetmedium, responsive to the solving.
 13. The system of claim 12 whereinenergy sources and energy detectors are configured in correspondingsource-detector pairs for the at least one scattering target medium andthe reference medium, and the modified perturbation equation has thefollowing form:${\left( {\delta\; I_{r}} \right)_{i} = {\left\lbrack \frac{I_{i} - \left( I_{0} \right)_{i}}{\left( I_{0} \right)_{i}} \right\rbrack\left( I_{r} \right)_{i}}};\mspace{14mu}{and}$W_(r) ⋅ δ x = δ I_(r) where δx is a vector of the relative changesbetween the known property of the reference medium and the correspondingunknown property of the at least one scattering target medium; W_(r) isa weight matrix describing the influence that each of a plurality ofvolume elements of the reference medium has on energy emerging at apoint on the reference medium; I_(r) is the vector of reference data;(I_(r))_(i) is an element of I_(r) for an ith source-detector pair ofthe reference medium; I_(i) is an element of the first vector for an ithsource-detector pair of the at least one scattering target mediumcorresponding to the ith source-detector pair of the reference medium;I₀ is the second vector; (I₀)_(i) is an element of I₀ for the ithsource-detector pair of the at least one scattering target medium;(δI_(r))_(i) represents a change in the vector of reference data for theith source-detector pair of the reference medium; and δI_(r) representsa change in the vector of reference data for all source-detector pairsof the reference medium.
 14. The system of claim 12 wherein thenormalizing of the difference between the first and second vectorscomprises determining a logarithm of a ratio of the first vector to thesecond vector.
 15. The system of claim 12 wherein energy sources andenergy detectors are configured in corresponding source-detector pairsfor the at least one scattering target medium and the reference medium,and the modified perturbation equation has the following form:${\left( {\delta\; I^{\prime}} \right)_{i} = {\ln\frac{I_{i}}{\left( I_{0} \right)_{i}}}};$${\left( W_{r}^{\prime} \right)_{ij} = \frac{\left( W_{r} \right)_{ij}}{\left( I_{r} \right)_{i}}};$δ I^(′) = W_(r)^(′)δ x where δx is a vector of the relative changesbetween the known property of the reference medium and the correspondingunknown property of the at least one scattering target medium; W′_(r) isa weight matrix describing the influence that each of a plurality ofvolume elements of the reference medium has on energy emerging at apoint on the reference medium; (W′_(r))_(ij) is an element of W′_(r) foran ith source-detector pair of the reference medium and a jth volumeelement of the reference medium; I_(r) is the vector of reference data;(I_(r))_(i) is an element of I_(r) for the ith source-detector pair ofthe reference medium; I_(i) is an element of the first vector for an ithsource-detector pair of the at least one scattering target mediumcorresponding to the ith source-detector pair of the reference medium;I₀ is the second vector; (I₀)_(i) is an element of I₀ for the ithsource-detector pair of the at least one scattering target medium; andδI′ is a vector representing a relative difference between the first andsecond vectors that is mapped onto the reference medium.
 16. Acomputer-implemented method for imaging of the properties of at leastone scattering target medium, wherein energy is directed into the atleast one scattering target medium using a plurality of energy sources,and the energy as it emerges from the at least one scattering targetmedium is detected using a plurality of energy detectors, the methodcomprising: generating first and second vectors of measured dataresponsive to readings from the energy detectors; generating a referencevector from a reference medium that models the at least one scatteringtarget medium, based on energy directed into the reference medium by aplurality of energy sources that model the plurality of energy sourcesof the at least one scattering target medium, and based on energydetected from the reference medium by a plurality of energy detectorsthat model the plurality of energy detectors of the at least onescattering target medium; calculating a normalized difference betweenthe first and second vectors of measured data according to a ratio of:(a) a difference between the first and second vectors, to (b) the secondvector; calculating a difference vector based on a product of thenormalized difference and the reference vector; solving a perturbationequation using the difference vector to obtain a vector representing aproperty of the at least one scattering target medium; and generating animage of the at least one scattering target medium based on the vectorrepresenting the property of the at least one scattering target medium.17. A medium on which computer code is embodied, the computer code beingexecutable by a computer to perform a method for imaging of theproperties of at least one scattering target medium, wherein energy isdirected into the at least one scattering target medium using aplurality of energy sources, and the energy as it emerges from the atleast one scattering target medium is detected using a plurality ofenergy detectors, the method comprising: generating first and secondvectors of measured data responsive to readings from the energydetectors; generating a reference vector from a reference medium thatmodels the at least one scattering target medium, based on energydirected into the reference medium by a plurality of energy sources thatmodel the plurality of energy sources of the at least one scatteringtarget medium, and based on energy detected from the reference medium bya plurality of energy detectors that model the plurality of energydetectors of the at least one scattering target medium; calculating anormalized difference between the first and second vectors of measureddata according to a ratio of: (a) a difference between the first andsecond vectors, to (b) the second vector; calculating a differencevector based on a product of the normalized difference and the referencevector; solving a perturbation equation using the difference vector toobtain a vector representing a property of the at least one scatteringtarget medium; and generating an image of the at least one scatteringtarget medium based on the vector representing the property of the atleast one scattering target medium.